A micro introduction to Bayesian modeling with BRMS

Julin Maloof

Frequentist vs Bayesian Statistics

Frequentist Statistics

  • What is the probability of the data given a hypothesis
  • Theory relies on data that is not observed
  • Accept or reject null hypothesis
  • Does not require use of priors

Bayesian

  • What is the probability of the hypothesis given the data
  • Does not rely on unobserved data
  • Evaluate ratio of probabilities for different hypotheses
  • Requires specification of prior beliefs (not as bad as it sounds!)

Prior probabilities

  • Bayesian methods require that we specify prior probabilities on the model parameters
  • This might sound like we are cheating
  • Instead we use relatively uninformative priors
  • As long as we have set relatively weak (uninformative) priors, the data will overwhelm the priors
  • (Although if we do have prior information we can use informative priors)

Sampling

  • Bayesian model fitting requires sampling different parameter values
  • For each set of parameters test, calculate the probability of the model, given the data
    • Similar to what we discussed two weeks ago with the nls function
  • Difference is that for nls the goal was just to find the “BEST” set of parameters
  • For Bayesian we need to sample enough to define the probability “landscape”
  • Typically several thousand sets of parameter values are sampled

Posterior probabilities

  • One output of Bayesian model fitting is a set of posterior probabilities
  • These will give us the distribution of possible parameter values for each model parameter, given the data and the priors
  • We can do further tests/comparisons once we have these posterior probabilities

Diagnostics

  • It is critically important to evaluate whether the sampling worked well
  • Typically four sampling “chains” are run
  • Rhat compares variance between and within chains. Should be < 1.05
  • ESS (Effective sampling size) should be > 100 per chain

Let’s try it

library(tidyverse)
library(brms)
dat <- read_csv("../output/height_data_clean.csv")

Please create two objects

  • datsmall that only retains observations from “2023-01-27”
  • datbcd that only retains BH, CC, and DPR observations from “2023-01-27”

Check the number of “cores” on our computer

  • Most computers have multiple “cores” or CPUs that can be used in parallel.
  • Let’s check what we have. We will use this on the next slide.
parallel::detectCores()
[1] 8

What are the default priors?

  • We will use what should now be a familiar formula to specify our model
  • We can ask what priors would be used by default for our formula and data
get_prior(height_cm ~ pop, 
          data = datbcd)
                  prior     class   coef group resp dpar nlpar lb ub
                 (flat)         b                                   
                 (flat)         b  popCC                            
                 (flat)         b popDPR                            
 student_t(3, 3.5, 2.5) Intercept                                   
   student_t(3, 0, 2.5)     sigma                               0   
       source
      default
 (vectorized)
 (vectorized)
      default
      default
  • They “b” class are the coefficients on fixed effects.
  • They have a flat prior (all values from -Infinity to +Infinity are equally likely)
  • We can do better than this while still not biasing the model (next slide)
  • Student_t is like a fat-tailed normal distribution.
    • The intercept is centered at 3.5,
    • sigma (standard deviation) is centered at 0

Student vs normal

Fit a model with height being predicted by pop

m1 <- brm(height_cm ~ pop,
          data = datbcd,
          prior = set_prior("normal(0,10)", class = "b"),
          sample_prior = TRUE,
          threads = 4) # set to 1 if you only have 1 core, or 2 if you have 2.
                       # don't go higher than 4, even if you have more than 4 cores.
  • We use an uninformative prior for “b”, the coefficient for differences between pops
    • We are setting this prior to be normally distributed, centered on 0, with a standard deviation of 10
    • Note “b” refers to any fixed-effect coefficients. If we had multiple fixed effects we could give them all the same prior, or we could specify them individually.
  • By default brms uses 4 chains threads = 4 means that we want each chain to run in parallel in a different CPU in our computer.
  • By default, each chain has 1,000 warm-up samples and 1,000 real samples

Fit a model with height being predicted by pop

m1 <- brm(height_cm ~ pop,
          data = datbcd,
          prior = set_prior("normal(0,10)", class = "b"),
          sample_prior = TRUE,
          threads = 4)
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 13.1.6 (clang-1316.0.21.2.5)’
using SDK: ‘MacOSX12.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   "-I/opt/R/arm64/include -I/opt/homebrew/include"    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000171 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.71 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.06 seconds (Warm-up)
Chain 1:                0.051 seconds (Sampling)
Chain 1:                0.111 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 6e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.055 seconds (Warm-up)
Chain 2:                0.053 seconds (Sampling)
Chain 2:                0.108 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 9e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.071 seconds (Warm-up)
Chain 3:                0.06 seconds (Sampling)
Chain 3:                0.131 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 1.3e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.091 seconds (Warm-up)
Chain 4:                0.096 seconds (Sampling)
Chain 4:                0.187 seconds (Total)
Chain 4:

The output shows us the progress of each chain

Model summary

summary(m1)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height_cm ~ pop 
   Data: datbcd (Number of observations: 155) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     3.22      0.12     2.98     3.47 1.00     3679     2887
popCC         1.45      0.25     0.96     1.93 1.00     4365     3309
popDPR        1.89      0.31     1.29     2.50 1.00     4094     3185

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.25      0.07     1.12     1.41 1.00     4271     2666

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
  • Estimate is the estimated height
    • Intercept is for the reference population (BH in this case)
    • popCC is the estimated difference between CC and BH
  • 95% confidence intervals show our confidence in the estimate.
    If they do not cross 0, then we are at least 95% confident that our differences are real
  • Be sure to check Rhat and the two ESS stats. What did we want these to be?

Model plots

plot(m1)
  • On the left we see the distribution of the posterior probability
    • Check for unusual pattern, bimodality, etc.
    • If the posterior is near the edge of the range of the prior, consider adjusting the prior
  • On the right we see the value for each parameter for each sample of each chain
    • Should look relatively uniform from left to right and chains should look similar

what are the other priors?

Our posterior plots had two coefficients that we did not specify priors for. What was used?

                  prior     class   coef group resp dpar nlpar lb ub
           normal(0,10)         b                                   
           normal(0,10)         b  popCC                            
           normal(0,10)         b popDPR                            
 student_t(3, 3.5, 2.5) Intercept                                   
   student_t(3, 0, 2.5)     sigma                               0   
       source
         user
 (vectorized)
 (vectorized)
      default
      default

Bayesian hypothesis testing

  • Is popCC really different from popBH?
  • Looking for
    • Post_prob to be near 1 (probability that hypothesis is true)
    • Evid.Ratio to > 3 (>3 = “moderate evidence”, > 10 = “strong evidence”)
hypothesis(m1, hypothesis = "popCC > 0")
Hypothesis Tests for class b:
   Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (popCC) > 0     1.45      0.25     1.04     1.86        Inf         1    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Getting posterior probability for an estimate

Sometime I prefer to do this by “hand”

posterior <- m1 %>% as_draws_df()
dim(posterior)
[1] 4000   12
head(posterior)
# A draws_df: 6 iterations, 1 chains, and 9 variables
  b_Intercept b_popCC b_popDPR sigma prior_Intercept prior_b prior_sigma lprior
1         3.2    0.74      1.9   1.3            3.42    16.5       3.920   -9.8
2         3.4    1.34      1.4   1.2           -0.48   -10.0       4.444   -9.8
3         3.1    1.79      1.8   1.3            3.09    16.4       3.284   -9.8
4         3.3    1.10      1.9   1.2            2.80    -4.8       1.350   -9.8
5         3.2    1.50      2.1   1.3           -1.50    -9.1       3.076   -9.8
6         3.2    1.44      2.0   1.4            3.85     2.8       0.077   -9.8
# ... with 1 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
# A tibble: 1 × 1
  CC_greater_than_zero
                 <dbl>
1                    1
  • All 4000 posterior samples estimate the the CC coefficient is > 0
  • This is very strong evidence the CC height is greater than BH height

Bayesian hypothesis testing

  • Is popCC the same height as popDPR?
  • Looking for
    • Post_prob to be near 1 (probability that hypothesis is true)
    • Evid.Ratio to > 3 (>3 = “moderate evidence”, > 10 = “strong evidence”)
hypothesis(m1, hypothesis = "popCC = popDPR")
Hypothesis Tests for class b:
            Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
1 (popCC)-(popDPR) = 0    -0.44      0.36    -1.15     0.28      19.59
  Post.Prob Star
1      0.95     
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

Fit a model with random effects:

  • Really pop and block should probably both be random, but we will keep pop fixed for now
m2 <- brm(height_cm ~ pop + (1|block),
          data = datbcd,
          prior = set_prior("normal(0,10)", class = "b"),
          sample_prior = TRUE,
          threads = 4)
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 13.1.6 (clang-1316.0.21.2.5)’
using SDK: ‘MacOSX12.3.sdk’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   "-I/opt/R/arm64/include -I/opt/homebrew/include"    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000185 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.85 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.34 seconds (Warm-up)
Chain 1:                0.352 seconds (Sampling)
Chain 1:                0.692 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 6e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.43 seconds (Warm-up)
Chain 2:                0.364 seconds (Sampling)
Chain 2:                0.794 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 2.7e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.537 seconds (Warm-up)
Chain 3:                0.291 seconds (Sampling)
Chain 3:                0.828 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 2.5e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.482 seconds (Warm-up)
Chain 4:                0.39 seconds (Sampling)
Chain 4:                0.872 seconds (Total)
Chain 4:

what priors were used for random effects?

prior_summary(m2)
                  prior     class      coef group resp dpar nlpar lb ub
           normal(0,10)         b                                      
           normal(0,10)         b     popCC                            
           normal(0,10)         b    popDPR                            
 student_t(3, 3.5, 2.5) Intercept                                      
   student_t(3, 0, 2.5)        sd                                  0   
   student_t(3, 0, 2.5)        sd           block                  0   
   student_t(3, 0, 2.5)        sd Intercept block                  0   
   student_t(3, 0, 2.5)     sigma                                  0   
       source
         user
 (vectorized)
 (vectorized)
      default
      default
 (vectorized)
 (vectorized)
      default

Check it (do your own)

summary

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height_cm ~ pop + (1 | block) 
   Data: datbcd (Number of observations: 155) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~block (Number of levels: 10) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.50      0.20     0.19     0.98 1.00     1212     1306

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     3.23      0.21     2.83     3.64 1.00     1567     1714
popCC         1.41      0.24     0.95     1.90 1.00     4199     2709
popDPR        1.85      0.28     1.28     2.41 1.00     4213     2675

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.18      0.07     1.05     1.33 1.00     3875     2936

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plots

Is it better than the model without the random effect?

m1 <- add_criterion(m1, "loo")
m2 <- add_criterion(m2, "loo")
loo_compare(m1, m2)
   elpd_diff se_diff
m2  0.0       0.0   
m1 -5.4       3.3
  • Preferred model is listed first
  • elpd_diff is how much worse a model is relative to the preferred model
  • Want this to be greater (ideally 2X greater) than se_diff to favor the more complex model